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Excitons are electron-hole pairs appearing below the band gap in insulators and semiconductors. 
They are vital to photovoltaics, but are hard to obtain with time-dependent density-functional 
theory (TDDFT), since most standard exchange-correlation (xc) functionals lack the proper long- 
range behavior. Furthermore, optical spectra of bulk solids calculated with TDDFT often lack the 
required resolution to distinguish discrete, weakly bound excitons from the continuum. We adapt 
the Casida equation formalism for molecular excitations to periodic solids, which allows us to obtain 
exciton binding energies directly. We calculate exciton binding energies for both small- and large- 
gap semiconductors and insulators, study the recently proposed bootstrap xc kernel [S. Sharma et 
al, Phys. Rev. Lett. 107, 186401 (2011)], and extend the formalism to triplet excitons. 
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I. INTRODUCTION 

Excitons arise from electron-hole attraction in gapped 
periodic systems such as bulk insulators and semiconduc- 
tors, as well as in many types of nanoscale systems, poly- 
mers and biomolecules.^'^ Bound excitons appear in op- 
tical spectra of extended systems as discrete absorption 
peaks below the quasiparticle gap, while continuum ex- 
citons enhance the band-edge absorption.'^ Excitons play 
an important role in photovoltaics, where photo-excited 
excitons propagate to heterojunctions and dissociate to 
yield currents. Although the phcnomcnological Wannicr 
model'^"^ describes excitons qualitatively well, it is not 
quantitatively suitable to be used in real applications 
where ab initio computation is required. 

The most important characteristic of bound excitons 
is their binding energy, defined as the difference between 
the quasiparticle gap and the excitation frequency of the 
exciton. The Bethe-Salpeter equation (BSE), a many- 
body method, is the standard way of calculating exciton 
binding energies in periodic systems,^ due to its accuracy. 
However, the scaling of the computational cost for BSE 
versus system size is not favorable, and the use of the 
BSE has therefore been limited to moderate system sizes, 
despite recent progress. 

Thanks to the balance of accuracy and computa- 
tional cost, density-functional theory (DFT) and time- 
dependent density-functional theory (TDDFT) are pop- 
ular ab initio methods for electronic structure and 
dynamics. ^"'^'"'^^ Instead of approaching the many-body 
problem directly, density-functional methods construct 
a noninteracting Kohn-Sham system with the same elec- 
tronic density as the original interacting system, which is 
much easier to solve than the original many-body prob- 
lem. Despite some additional difficulties for periodic sys- 
tems (in particular, the severely underestimated gap), 

TDDFT methods are gaining popularity in solid-state 
physics.6'13,14 

TDDFT is a formally exact theory for electron dynam- 
ics, but in practice the exchange-correlation (xc) kernel 



must be approximated. It has been notoriously diffi- 
cult to get excitons in TDDFT: ^'"'"^^ local and semilo- 
cal xc kernels that work well in finite systems do not 
yield bound excitons in solids, since they lack a long- 
range part.^" The recently proposed long-range correc- 
tion (LRC) xc kcrnel^°"^^ allows bound excitons to be 
obtained from TDDFT, but empirical input is required. 

Aside from the difficulty to find good xc kernels, there 
is another problem. TDDFT approaches for periodic sys- 
tems typically calculate the optical spectrum via the di- 
electric function; but exciton binding energies in semi- 
conductors are usually in the mcV range, which means 
that bound excitons require a high frequency resolution 
to be distinguished from the continuum. This makes the 
calculation numerically demanding. One could increase 
the frequency resolution near the region of interest, but 
this requires knowing the exciton binding energies before- 
hand. As a consequence, most existing TDDFT studies 
of excitons are either for materials with strongly bound 
excitons far away from the band edge such as LiF or Ar,^'' 
or describe the enhancement of the band-edge continuum 
spectrum due to excitonic effects. ^'^'^ 

We recently proposed an alternative TDDFT approach 
for obtaining excitonic binding energies directly. ^^^^^ 
The approach was applied to one-dimensional model 
systems, where we showed that TDDFT within the adi- 
abatic approximation can yield more than one exciton if 
local-field effects are included. We also considered sev- 
eral bulk solids and found that TDDFT, using xc kernels 
with appropriate long-range behavior, can yield excitonic 
binding energies in the right range. However, this ear- 
lier work remained somewhat inconclusive due to several 
simplifications (most notably, a two-band approximation, 
a rather small fc-space grid, and a real-space representa- 
tion of the xc kernel which resulted in a loss of accuracy) . 

In this paper we present a systematic computational 
study of the lowest excitonic binding energies in com- 
mon zincblende and wurtzite semiconductors as well as 
in large-gap insulators. We extend our earlier work^^ in 
several ways: we go beyond the two-band approxima- 
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tion and include, in principle, an arbitrary number of 
bands; this can be viewed as the solid-state analog of 
the Casida approach for molecular excitation energies. 
Furthermore, we extend the formalism to include triplet 
excitons, and we test the so-called bootstrap xc kernel.^® 
Atomic units (e = = me = l/(47reo) = 1) are used 
throughout this paper unless mentioned otherwise. 



II. THEORETICAL BACKGROUND 

In semiconductors, the binding between an electron 
and a hole is usually weak: such electron-hole pairs are 
designated as Wannier excitons. The Wannier model^"^ 
describes such excitons in analogy to positronium, where 
the effect of the material environment is introduced by 
the effective mass and the dielectric constant. The Wan- 
nier equation for excitons is given by 

^,(r) =£;,^,(r), (1) 
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where = (jti^^— m^T^)^^ is the reduced effective mass, 
e is the dielectric constant of the material, and E^, and 
"tpi, are excitonic binding energy and wave function, re- 
spectively. The binding energy is the most important 
property for excitons, defined as the difference between 
the quasiparticle gap and the excitonic excitation energy. 
Despite the simplicity of Eq. (1), its exciton binding 
energies can be fairly accurate for common semiconduc- 
tors such as GaAs^^ and Cu20,^° since the model can 
be derived as an approximation to the BSE many-body 
theory'^^ (assuming the effective Bohr radii of excitons 
are much greater than the lattice constant). 

Excitons in large-gap materials (such as LiF and Ar) 
arc strongly bound and localized within a single crystal 
unit. The Bohr radii of these so-called Frenkel excitons 
are small, so the Wannier model does not describe Frenkel 
excitons well. However, from the point of view of an ab 
initio electronic structure theory, there is no conceptual 
difference between Wannier and Frenkel excitons, as they 
all are just excitations of the many-body system. The col- 
lective character of the excitons distinguishes them from 
other excitations, i.e., they arise from superpositions of 
many single-particle excitations. 

Aside from their collective quasiparticle character, ex- 
citons are normal optical excitations. TDDFT has been 
successful in treating excitations in finite systems, and its 
use for periodic systems is increasing. TDDFT solves a 
non-interacting time-dependent system described by the 
time-dependent Kohn-Sham equation: 



.^0(r,O = 



V2 

2 



-^^ Wcxt(r,t) + VH(r,<) 



+ Wxc(r,i) 



(2) 



where ?;oxt and wh arc the external potential and the 
Hartree potential, respectively, Wxc is the exchange- 
correlation (xc) potential, and is a time-dependent 



Kohn-Sham orbital, fxc is defined as the one-body mul- 
tiplicative potential with which the solution of Eq. (2) 
reproduces the density of the interacting system, and it is 
the only part that needs to be approximated in practice. 
One can obtain information about the excitations in the 
interacting system by propagating Eq. (2) under an ex- 
ternal perturbativc potential, '^^ but the more convenient 
approach is to work in the frequency domain from the 
beginning. 

The linear response function'^'^ x = ^n/Sv^j^t describes 
the first-order density change caused by a change in the 
external potential, and thus determines the optical spec- 
trum. In TDDFT, the response function in reciprocal 
space is obtained, in principle exactly, as 



XGG'(q,w) = 



G' 



J2^'^-GiG3i<l,Uj) 



Ga 



/] 



Hxc.GaG 



Xs^G"G'{^,(^), (3) 



GG" 



where Xs is the linear response function of the Kohn- 
Sham system, "'^^ /hxc = /h + /xc with the Hartree kernel 
/h = Svn/Sn (in reciprocal space, /h = AttSgg' /\(l+G\), 
and the xc kernel /xc = Sv^c/^n-. All quantities in Eq. 
(3) are matrices in momentum space, where q belongs to 
the first Brillouin zone, and the G's are reciprocal lattice 
vectors. 

The optical absorption spectrum of a periodic system 
is described by the macroscopic dielectric function em-^ 
One calculates em from x by 



eM(w) = lim 



1 



q^o 1 + 47rxoo(q,w)/g2 



(4) 



Calculating em on a frequency grid yields the spectrum. 

The so-called head (G = G' = 0) of the xc kernel 
gives the largest contribution to the change from Xs to 
X, and the contributions from bigger G's decay rapidly. 
Thus the sums in Eq. (3) can be restricted to a small 
number of reciprocal lattice vectors, which reduces the 
computational effort significantly. 

One can select the frequency uj in Eq. (3), so the cal- 
culation can be focused on the region of interest instead 
of the entire spectrum. It might appear that in Eq. (3) 
one can choose any frequency resolution, but it is im- 
plicitly limited by the number of Kohn-Sham excitations 
included in x^. A part of the spectrum is directly ob- 
tained from this approach, but there is no way of know- 
ing which Kohn-Sham excitations contribute to a specific 
peak. Considering the continuum nature of the spec- 
tra of periodic systems, these details are of course rarely 
needed. For excitons, however, the binding energies are 
not explicitly given in this approach, unlike in the much 
simpler Wannier model. Due to the discrete nature of 
bound excitons, the Kohn-Sham excitation composition 
is useful for interpretations; but this information is not 
available in this approach. 

In finite systems the low-lying excitations are discrete, 
and there is a more efficient way to obtain them than 
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scanning the frequency range with Eq. (3). The idea 
is to describe electronic excitations as eigenmodes of 
the system. '^^ The Casida equation^'^ then transforms 
the TDDFT Hnear-response equation into the transition 
space spanned by single-particle Kohn-Sham excitations; 



A B 
B A 



-1 




(5) 



(6) 



where the matrix elements of A and B are 

with the Hartree-exchange-correlation (Hxc) matrix F^xc 
for spin-unpolarized systems 



F 



{ij)(mn] 



d\ / d\' cb,{r)cf,*{r)Uxc{r,r',Lj) 



x<p*m{r')Mr')- (7) 

The factor of 2 in Fhxc accounts for the spin and the 0's 
are the ground-state Kohn-Sham orbitals. The indices i, 
j, m, n all represent full sets of quantum numbers, where 
i, m denote occupied orbitals and j, n denote unoccupied 
orbitals. 

Most of the currently available xc kernels are frequency 
independent, in which case Eq. (5) becomes an eigen- 
value problem. The explicit matrix formulation of Eq. 
(5) is suitable for discrete excitations in finite systems. 
The excitation frequencies of the system are explicitly 
given by the eigenvalues uj. The eigenvector X together 
with Y describes how the Kohn-Sham excitations com- 
bine to form the excitation in the real system. The op- 
tical spectrum can be calculated with X and Y.^ The 
widely used Tamm-Dancoff approximation (TDA) sets 
the matrix B to zero and hence neglects the correlation 
between excitations and de-excitations. We have shown 
that the TDA can often be a better choice for excitons 
than the exact calculation with Eq. (5),^^ and we there- 
fore employ the TDA throughout this paper. 

The difficulty of obtaining excitons in TDDFT mainly 
comes from the requirement on the xc kernel for peri- 
odic systems, that it needs a behavior in the q — > 
limit. This long-range behavior is necessary in order to 
produce non-zero head (G = G' = 0) and wing (G = or 
G' = 0) contribution to x in Eq. (3) and to Fxc in Eq. 
(7). In periodic systems these contributions dominate 
over the so-called local-field effects [contributions from 
the body (G 7^ and G' 7^ 0) of /xc]- Common local and 
semi-local xc kernels, such as the adiabatic local-density 
approximation (ALDA), lack this long-range behavior. 
Several recently proposed functionals (such as the empir- 
ical long-range correction^^, the non-empirical bootstrap 
kerneP* and the non-empirical meta-GGA kernel^^) have 
the correct long-range behavior; hence, they are promis- 
ing choices for the accurate calculation of excitons in 
TDDFT. Due to its discrete nature, the Casida equation 
approach is rarely used for periodic systems. For calcu- 
lating cxciton binding energies, however, it is preferable 
over the usual response-function approach. 



III. METHOD 

In this section we present the details of how we ap- 
ply our TDDFT approach for calculating exciton bind- 
ing energies. Within the TDA and using the adiabatic 
approximation for the xc kernel, Eq. (5) becomes 



(ij)(mn) 

HXC 



{ran) 



^SimSjniej - ei) + Fi 

(mn) 

(8) 

We only consider optical excitations which have no mo- 
mentum transfer. Thus only the q = part of the xc 
kernel is involved. The xc matrix element in reciprocal 
space for spin-unpolarized systems is given by 



j-i(ijk)(mnk') 
^xc 



/xc,GG'(q = 0) 



JG-r\ 



k) (mk' 



-iG 'r 



nk'),(9) 



where V is the volume of the crystal, k, k' are Bloch 
wavevectors of orbitals, and i, j, m, n are band in- 
dices. Since the matrix element (jk |e'*^ ""| zk) vanishes 
as G — >■ 0, the contributions of the head (G = G' = 0) 
and of the wings (G = or G' = 0) need to be evaluated 
analytically, /xc must diverge as for the head and as 
q^^ for the wings for these to have a nonzero contribu- 
tion. In this work, we use the ABINIT pseudopotential 
band structure code for the Kohn-Sham ground state. 
Due to the presence of pscudopotentials, the above ma- 
trix element for G = must be replaced by'^'' 



iG-i 



— )• 



{jk\p^i[f,V^,]\ik) 

ejk — Elk 



(10) 



where p is the momentum operator, f is the position oper- 
ator, and is the nonlocal part of the pseudopotential. 



A. XC kernels 

Wc use the following adiabatic xc kernels in our study: 
the long-range correction, the bootstrap kernel, and 
the PGG kernel for singlet excitons, '^^''^^ and the hybrid 
kernel by Burke et al. for triplet excitons. 

The long-range correction kernel is defined as 



/x^|^^(q,G,G') 



a 



|q+G| 



r^G', (11) 



where a is a material-dependent parameter. In Rcf. 21, 
an empirical formula was proposed for determining a: 



a = 4.615e;^^ - 0.213, 



(12) 



where £00 is the high-frequency dielectric constant. The 
purpose of this empirical formula is to reproduce the con- 
tinuum spectrum; here, we test its effect on the exciton 
binding energy, which would be too small to be resolved 



4 



in common response-function calculations.^^ For compar- 
ison, we will also fit a with respect to the experimental 
exciton binding energies. 

The PGG kernel is a real-space kernel approximating 
the exact exchange kernel. In real space, it is defined 



as' 



38,39 



/^"(r,r') = -- 



|r — r'l 7i(r)n(r') 



(13) 



where n is the ground-state electronic density. We con- 
vert the PGG kernel into reciprocal space for its use in 
Eq. (9). The Kohn-Sham orbitals in Eq. (13) have the 
Bloch form. 



1 



^^.k(r)e"'•^ 



(14) 



where A^ceii is the number of unit cells in the crystal, and 
Wik(r) is the Bloch function, /xc?*^ can then be written 
as 



occ. occ. 



Op-i(k-k')-(r-r') 

/--(r, r') = -l^^ — if,.,„k' (r, r'), 

(15) 

where -ff^kink' (r, r') is periodic within one unit cell and 
defined as 



H. 



■ikmk 



(r,r') 



(16) 



The Fourier transform of f^^'^ yields 



xc 
occ. occ. 



/--(q,G,G') = 4EEE 



Stt 



X-f^'ikmk'(G — Go, G' - Go), (17) 

where H is obtained by numerical Fourier transform of 
expression (16) within one unit cell. For simplicity, we 
ignore the local-field effects and only use the head of the 
PGG kernel, which is given by 



8^ ^ i/,kmk(0,0) 



(18) 



The PGG kernel is orbital dependent and involves a sum 
over all occupied orbitals. It is not obvious whether the 
pseudopotential formalism is directly compatible with 
this kernel since it does not properly include the core 
states; this will be the subject of future study. For the 
time being, we use all the occupied pseudo-bands for con- 
structing the PGG kernel. 

The so-called bootstrap kcrncl^^ is a recently proposed 
non-empirical adiabatic xc kernel, designed to be able to 
treat excitons: 



/x™"''(q,G,G') 



ts,sym 



(q,G = G' = 0) 



, (19) 



where /xc,sym, Csyln Xs.sym are the xc kernel, inverse 
dielectric function and Kohn-Sham linear response func- 
tion in their symmetric forms, respectively, defined as 

/xc..ym(q, G, G') = I'G ^'(q)/xc(q, G, G')i'G'^'(q!?0) 



e,-! (q, G, G') = t>^i/^(q)e-i(q, G, G')«^^,^(q),(21) 



Xsym(q, G, G') = 4^'(q)x(q, G, G')t;G'''(q), (22) 

where WG(q) = 47r/|q + G|^ is the Coulomb potential. 
e~\^ in Eq. (19) is calculated as 



e"^ = 1 -I- y 

^sym ' A 



sym 



= 1 + (1 - Xs,sym/i 



HXC,sym 



(23) 



s,symi 



where all quantities in Eq. (23) are matrices in G and 
G'. Equations (23) and (19) are iteratively evaluated 
until self-consistency is achieved. Since the head and the 
wings of /xc diverge as q — > (which is important for 
excitons), /xc cannot be used directly in the iteration due 
to numerical difficulties. The symmetric forms ensure 
that no troubling singularities are involved. 

So far, we have neglected the spin. In principle, a 
noncoUinear spin formulation for the xc kernel^^'^^ is 
needed to treat singlet and triplet excitations on the same 
ground. For spin-unpolarized systems, however, one can 
define singlet and triplet xc kernels as 



/: 



singlet 

xc 



tt 
xc 



fU 
./xc 



/: 



triplet 

xc 



ftt 

Jxc 



/: 



u 

xc 



(24) 



where /xc = Sv^ca/Sna-i is the spin-dependent xc kernel, 
eft 



Only /xi and /xc are involved because /xc = ./xc and 



/xc = /xc in spin-unpolarized systems. One calculates 
singlet and triplet excitations by performing two separate 
TDDFT calculations with /xc^'°* and /xcP'''*(the triplet 
calculation does not include the Hartree kernel). In the 
spin-dependent case, the important property for excitons 
is the exchange splitting Ax, defined as 



Ax = K 



triplet 



-e: 



isinglct 



(25) 



The bootstrap xc kernel was originally developed for 
the spin-independent case. Since it is not defined 
via a functional derivative, there is no unique way to 
make it spin dependent. We propose a plausible spin- 
dependent generalization of the bootstrap kernel in such 
a way that the singlet result does not change during 
the self-consistent procedure. The idea is to replace 



X's 



y^{q, G = G' = 0) in Eq. (19) by a matrix Ma-a' 

fxc°sym^aa' G, G') = ^ <^sym,aa" i^' G, G") 



a"G" 



xAV,,(q,G",G'), (26) 

where the spin-dependent inverse dielectric function e^^ 
for a spin-unpolarized system can be defined as 



'Xo 



(27) 
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The Coulomb interaction is spin-independent {va 
so in Eq. (27) becomes 



sym.(T(T 



= s„ 



vsym 



V-sym 



(28) 



To ensure that fxc^^°^ 



/bootstrap i • j rj- ir 

xc during and atter selt- 



consistency, it is straightforward to show that the matrix 
Maa' must satisfy the foUowing conditions: 

M,^,,(q, G = G') = -Af,=,,(q, G = G') 

+ 2x.:lym(q,G-G' = 0),(29) 
M<,^,,(q, G ^ G') = -A-W(q, G ^ G')- (30) 

This ensures that the singlet xc kernel reproduces the 
bootstrap kernel, but it also has consequences for the 
triplet xc kernel. If the M^^i matrix is the same for all 
steps of the iteration, it is easy to show that during the 



iteration /; 



triplet 



would not change, and it can be deter- 



mined without any iterative procedure as 

/xc'!'sym(q, G, G') = i [il/n(q, G, G') - Mniq, G, G')] . 

(31) 

jtiipiet determined self-consistcntly only if Mo-o-' de- 
pends on /xc of the previous iteration step. 

Burke, Petersilka, and Gross^° proposed a hybrid spin- 
dependent xc kernel, where 



/■hybrid _ rPGG 
/xc,tt /xc,tt' 

j^hybrid j-ALDA 

J xc,T4- J xc,T4- ■ 



(32) 
(33) 



This hybrid kernel yields rather accurate singlet-triplet 
splittings as well as excitations frequencies in finite sys- 
tem. Inspired by this, we propose to calculate the singlet- 
triplet splitting of excitons with a similar hybrid kernel. 
Since the bootstrap kernel is a parameter-free xc kernel 
that has good performance for exciton binding energies 
in several tested systems, we use the bootstrap kernel as 
the singlet xc kernel: 



/^bootstrap n ^bootstrap 

•/xCjI-t ^Jxc 



/•ALDA 

JXCjtJ. ' 



xc,"t"4- 



xc,ti 



(34) 
(35) 



The singlet exciton binding energy retains the result of 
the bootstrap kernel. This modified hybrid kernel is ef- 
fectively an instance of the generalized spin-dependent 
bootstrap kernel, where 

i\,r _ i-bootstrap yALDA , 0^—1 ('ifi'] 

and Mfi is determined by Eqs. (29) and (30). 



B. Computational aspects and numerical details 

Let us now briefly summarize the computational as- 
pects of our exciton calculations; more technical details 
will be presented elsewhere. 



An LDA ground-state calculation is carried out with 
the ABINIT pseudopotential code.'^^ The ground-state 
band structure and Bloch functions arc taken as the in- 
put to our TDDFT calculations. Due to the absence of 
the derivative discontinuity, the LDA Kohn-Sham gap is 
too small to approximate the quasiparticle gap.^^ This 
introduces big errors in the corresponding TDDFT cal- 
culation, since adiabatic xc kernels cannot change the 
gap.^^ A frequency-dependent kernel together with a 
good ground-state xc functional would be required to 
fully treat the gap problem in TDDFT. However, since 
our focus is on the exciton binding energies and not on 
the gap itself, we shift the gap to its corresponding ex- 
perimental value by applying a simple scissor operator^^ 
to the conduction bands, and we apply a corresponding 
correction to the momentum matrix elements.'*^ For xc 
kernels that explicitly depend on the density, we use the 
pseudodensity with the correction described in Ref. 47. 

The Fxc matrix in Eq. (9) is represented in the transi- 
tion space of Kohn-Sham excitations. It has the dimen- 
sion iVv X iVc X iVk, where v stands for valence bands, c 
stands for conduction bands, and iVk is the number of k- 
grid points in the Brillouin zone. To achieve convergence, 
the dimension must be large, so a Casida-equation-type 
calculation is computationally more demanding than the 
usual response- function calculation; in particular, storage 
of a large matrix is required. For this reason, Casida- 
equation-type calculations are not usually done for peri- 
odic systems. In contrast, a response-function calculation 
processes matrices in reciprocal space; the dimension of 
the matrix equals the number of reciprocal lattice vectors 
used in the calculation, which is much smaller than the 
size of the Fxc matrix [Eq. (9)] required for convergence. 

According to the Wannier model, excitons are dom- 
inated by single-particle (Kohn-Sham) excitations near 
the band edge for direct-gap solids. Although strictly 
speaking the Wannier model refers to quasiparticle bands 
instead of Kohn-Sham bands, the Kohn-Sham band 
structure is similar to the quasiparticle band structure 
near the Fermi level (aside from having the wrong gap).^^ 
Thus for direct-gap solids, only the highest valence bands 
and the lowest conduction bands are needed in the cal- 
culation (including degeneracies at the F point). 

With only a few bands used, the dimension of the ma- 
trices in Eq. (5) is vastly reduced, so the eigenvalue 
problem can be solved with acceptable cost. For the 
same reason, we also only include excitations between 
these bands in the Kohn-Sham response function used 
in the bootstrap kernel. Although this restriction to a 
few bands near the band gap is sufficient for describ- 
ing excitons, this would yield unsatisfactory continuum 
spectra; but those are not our concern in this work, since 
they can be calculated much better with the standard 
response-function approach of TDDFT. 

To make the numerical problem more manageable for 
resource-limited environments, the storage of the large 
matrices can be distributed among different nodes of a 
computer cluster, and we use the ScaLAPACK^'^ package 
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to manipulate the distributed matrix. Since we only need 
the first exciton binding energy, we can further simplify 
the eigenvalue problem by using an iterative eigenvalue 
solver such as FEAST^^ to converge only the eigenvalue 
for the exciton instead of all the excitations. Iterative 
eigenvalue solvers require a predefined range of the de- 
sired eigenvalues, which would require knowledge about 
the result before it is obtained. But in our case the range 
of frequency is conveniently chosen from zero to the band 
gap, since excitonic excitations are always below the band 
gap. With all the techniques described in this section, 
Casida-equation-like calculations for periodic solids be- 
come practically feasible. 

We check convergence against the following parame- 
ters: the number of bands in the calculation, the k-grid 
for the ground state, and the reciprocal lattice vectors in 
constructing the matrix in Eq. (8). 

For all systems under consideration, the highest va- 
lence band has p-character and is triply degenerate at 
the F point, while there is no degeneracy at the F point 
for the lowest conduction band. Thus we only include 
3 valence bands and 1 conduction band in Eq. (8) and 
in the Kohn-Sham response function Xs for the boot- 
strap kernel. We tested the convergence with respect to 
number of bands, and we found that the effect in the 
exciton binding energy from increasing the number of 
valence bands is much smaller than that of increasing 
the number of G- vectors. We test convergence against 
the TDDFT response-function calculations in Ref. 28, 
which by nature include much more bands than we use 
here. The excitons are strong enough in Ar and LiF to 
be resolved with the response-function approach, and our 
calculation with three valence and one conduction bands 
produces exciton binding energies very close to those re- 
ported in Ref. 28 (see Table I). We thus conclude that 
our few-band approach is sufficient for exciton binding 
energies. 

We find that for strongly bound Frenkel excitons in in- 
sulators, a small k-grid such as a 10 x 10 x 10 Monkhorst- 
Pack grid^° is sufficient for the convergence of the binding 
energy. In the case of weakly bound Wannier excitons in 
semiconductors, the convergence is much slower and re- 
quires at least a 18 x 18 x 18 grid. This result is not 
surprising, since Frenkel excitons are local in real space, 
which means they are diffuse in reciprocal space, and 
therefore require less resolution to be well-described in 
the reciprocal space than Wannier excitons. 

Several previous TDDFT studies (within the response- 
function approach) state that convergence is achieved for 
semiconductors with grids like 15 x 15 x 15.^^'^^ How- 
ever, in these cases convergence is to be understood with 
respect to the continuum spectra instead of the exciton 
binding energies. We find that with a 15 x 15 x 15 grid, the 
maximum relative error in exciton binding energies for 
semiconductors studied in this work is 139%. We there- 
fore calculate zincblende materials with a 18 x 18 x 18 
grid, wurtzite materials with a 20 x 20 x 20 grid, and 
insulators with a 10 x 10 x 10 grid. The ground-state cal- 



culation with a bigger k-grid is feasible since the number 
of effective k-points can be greatly reduced by symmetry, 
but this is not the case for TDDFT, where the k-points 
in the entire first Brillouin zone must be included. 

Excitons are usually described as electron-hole pairs 
and are modeled with parabolic bands. This may suggest 
that one can get away with using only k-points near the 
F-point (k = 0), which would greatly decrease the size of 
the problem. This approach was employed by Rohlfing 
and Louie in a BSE study.^^ In a TDDFT context, how- 
ever, for GaAs including 62% of the k-points (centered 
at the F-point) induces a 74% relative error, and even in- 
cluding 95% of the k-points still induces a 5% error, and 
the benefit of calculation speed is diminishing. We con- 
firm the result with BSE and TDDFT calculations for an 
ID Kronig-Penney model, where we find that only in- 
cluding fc-points near fc = induces a relative error three 
times larger in TDDFT than in BSE. This can be under- 
stood by analyzing the coupling matrix of TDDFT [Fxc 
in Eq. (9)] and BSE (Ref. 26): the equivalent object of 
Fxc in BSE is dominated by its diagonal (k = k') part, 
so only taking the k-points near the F-point still retains 
its shape. In TDDFT this is not the case. Excitons are 
known to have collective character, but this discrepancy 
between model systems and excitons in real materials 
points to a surprisingly high degree of collectiveness when 
represented with single-particle Kohn-Sham excitations, 
since Kohn-Sham excitations contribute over the entire 
Brillouin zone. 

The head of the xc kernel makes the largest contribu- 
tion in Eq. (8); the other contributions can be usually 
ignored without loss of accuracy. We find that this is 
also true for the exciton binding energy. We checked the 
error in the exciton binding energy introduced by only 
using the head of the xc kernel, and the error is about 
1% in most cases, and less than 5% throughout. Thus 
we only use the head in most calculations, but in cases 
where the head vanishes (such as the ALDA part included 
in the hybrid kernel), we include G- vectors with length 
up to 2Go, with Go being the longest reciprocal cell vec- 
tor. The construction of the bootstrap kernel involves 
matrix operations in the reciprocal space, so we also use 
G-vectors up to 2Go for it. After the bootstrap kernel 
is calculated, we only use its head in Eq. (8). It should 
be noticed that although this still yields acceptable ac- 
curacies for the exciton binding energy, only using the 
head would not produce more than one exciton. If one 
needs an excitonic Rydberg series, the wings and body 
contribution of the xc kernel are necessary. 

For the bootstrap kernel there is an additional conver- 
gence issue related to the number of bands that need to 
be included in the Kohn-Sham response function Xs- For 
large-gap insulators, we only need to use four bands in 
the calculation of Xs', For zincblende semiconductors, we 
have to use a total of 60 bands to achieve convergence. 
For wurtzite semiconductors, we use 10 valence bands 
and 40 conduction bands in the calculation. 
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TABLE I. Lowest singlet exciton binding energies ^^'"sict singlet-triplet splittings Ax, calculated with TDDFT using 
various different xc kernels, and compared to experimental values from the literature. A star (*) means that no bound exciton 
was obtained from the calculation, "n.c." means that no calculation was performed. 



Material'^ 


GaAs 


P-GaN 


Q-GaN 


CdS 


CdSe 


Ar 


Ne 


LiF 


Cjxp. gap (e\ ) 


1 


o.o 






1 74 




91 ^il 


1 4 9n 


rjxp. _C/jj 


^ 97mpV 


OR nmpV 


Of) ztmpV 


9e OttipV 


1 rimpv 

-L .Ullic V 


1 qripV 


4 nspv 


1 fipV 




Q fil //pV 








49 78(/pV 


0.16eV 


0.14eV 




LRC empirical 


0.211 


0.6578 


0.6496 


0.6448 


0.4721 


2.697 




2.191 


LRC empirical e;^'"^'^* 


0.8580meV 


0.5143meV 


* 


0.5131meV 


1.405meV 


0.3043meV 




1.136meV 


LRC fit Q 


0.595 


2.409 


3.6285 


4.244 


2.144 


21.45 


96.5 


9.5 


/■bootstrap rrisinglct 

/xc 


0.3318meV 


0.1992meV 


* 


0.4610meV 


0.8947meV 


2.156eV 


6.225eV 


1.547eV 


Corresponding LRC a'^ 


0.08836 


0.3048 


0.2147 


0.5895 


0.3183 


22.6324 


126.673 


9.32326 


PGG Bf's'<=« 


* 


n.c. 


n.c. 


n.c. 


n.c. 


* 


n.c. 


* 


Corresponding LRC a"^ 


5.820 X lO"-"^ 


n.c. 


n.c. 


n.c. 


n.c. 


3.924 X 10"" 


n.c. 


3.851 X 10""* 


hybrid^ Ax 


37.65AteV 


40.5lAteV 


* 


n.c. 


n.c. 


0.03297eV 


0.01287eV 


0.5041meV 



Unless otherwise mentioned, zincblende materials are calculated with 18 X 18 X 18 Monkhorst-Pack k-point grid, wurtzite materials are 
calculated with 20 X 20 X 20 grid, and solid Ar, solid Ne, and LiF are calculated with 10 X 10 X 10 grid. 
Experimental data from Rcfs. 52-59. 

Experimental data from Refs. 57, 59-61. No experimental data available for GaN, CdS, and LiF. 

Calculated using Eq. (12). The data is not available for solid Ne. 

The head of the xc kernel has the same form as the LRC. 
' Local field effects included, 
s See Eqs. (34) and (35). 



IV. RESULTS 

We consider several common direct-gap zincblende 
(GaAs, /3-GaN) and wurtzite (a-GaN, CdS, CdSe) semi- 
conductors, as well as insulators (LiF, solid Ar, solid Ne). 
Results are presented in Table L The first three rows give 
experimental data on the band gap,^"*'^^"^^ the binding 
energy ^^'"sict ^£ ^^^^ lowest singlet exciton, ^^"^^ and the 
singlet-triplet exchange splitting Ax.^'''^^"®^ The remain- 
ing rows of Table I show the results of our calculations. 

The exciton binding energies calculated with the LRC 
kernel are generally found to be significantly too small if 
the empirical a is used, see Eq. (12); for a-GaN, there 
is not even a bound exciton. We therefore determine the 
LRC a parameter for the exciton binding energies by fit- 
ting to experimental data for ^™siet^ ^^^^ that 
they are quite different from the empirical formula for 
a. Ref. 21 argued that a must be proportional to e^, 
but we cannot confirm such a linear fit for a with our re- 
sults. The sensitivity of the exciton binding energy with 
respect to a varies a lot for different materials. Though 
the empirical formula for a [Eq. (12)] was originally not 
developed to give accurate exciton binding energies, we 
find that calculations with empirical a still yield bound 
excitons (except for a-GaN). This may explain, at least in 
part, why these parameters lead to quite accurate spectra 
in the vicinity of the gap. 

We next consider two nonempirical kernels, bootstrap 
and PGG. We find that the self-consistent procedure of 
the bootstrap kernel is very stable: even if we start the 
iteration with a different /xc (instead of starting with no 



/xc), we always converge to the same bootstrap kernel. 
We confirm that the bootstrap kernel produces accurate 
exciton binding energies for Frenkel excitons in Ar and 
LiF, as reported in ReL 28. For Ne the bootstrap kernel 
overbinds by about 50%, but it still yields the correct 
order of magnitude for ^'^s''^* , Por Wannier excitons in 
the studied semiconductors, however, the bootstrap ker- 
nel fails to differentiate between different materials and 
in all cases yields exciton binding energies that are too 
low. On the other hand, the excitonic enhancement of the 
continuum spectrum is reported to be well-described by 
this kernel. This is understandable since the correspond- 
ing LRC a parameters for semiconductors are close to 
those given by the empirical formula Eq. (12). 

The performance of the PGG kernel, which works well 
in finite systems, is disappointing: it does not produce 
any bound excitons at all, despite having a nonzero 
head contribution. The PGG kernel is an exchange- 
only kernel, and is known to bind quite strongly in finite 
systems, so it is surprising that it does not yield any 
bound exciton in the cases we tested. One possible reason 
is that the pseudopotential treatment is not compatible 
with the explicit orbital dependence in the PGG kernel, 
since the contribution from core orbitals cannot be sys- 
tematically included. Also it should be noted that while 
periodic systems are dominated by the head of the xc 
kernel, there is no corresponding effect in finite systems. 
This is because in finite systems the electron dynamics 
can be viewed as coming entirely from local-field effects. 
Thus, the strongly attractive nature of the PGG kernel 
in finite systems would at most translate into a strong 
body of the xc matrix in periodic systems (which, how- 
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ever, is irrelevant for excitons), but does not necessarily 
guarantee a strong head. This is indeed confirmed by 
calculating the LRC a that corresponds to the PGG ker- 
nel, which turns out to be orders of magnitude too weak 
(see the second-to-last row of Table I). 

As wc discussed in Section III A, not many long-range 
spin-dependent xc kernels for treating triplet excitons are 
known. Our hybrid kernel [Eqs. (34) and (35)] can be 
viewed as a special case of the generalized bootstrap ker- 
nel in Eq. (26), which performed well for singlet excitons. 
We find that the exchange splitting is of the right order 
of magnitude for GaAs, although somewhat too large. A 
similar Ax is found for /3-GaN, but there is no experi- 
mental data available for comparison. Since a-GaN does 
not have a bound exciton with the bootstrap kernel, it is 
not possible to obtain a well-defined Ax with the hybrid 
kernel. No calculation was performed for CdS and CdSc 
due to excessive memory requirements. Finally, for the 
strongly bound excitons in the large-gap insulators Ar, 
Ne, and LiF, Ax comes out significantly too small. 

Considering the fact that the singlet-triplet exchange 
splitting is essentially treated on an ALDA level, it is 
perhaps not surprising that the results for Ax with the 
hybrid kernel are not terribly accurate. Clearly, more 
sophisticated spin-dependent xc kernels for singlet and 
triplet excitons need to be developed. 

V. CONCLUSION 

In this paper we have introduced an alternative 
TDDFT approach for calculating excitonic binding en- 
ergies in solids. We present the first converged Casida- 
equation-type TDDFT calculations for several materials, 
showing that such calculations are feasible for real peri- 
odic bulk systems. The approach yields exciton binding 
energies directly, rather than the optical spectrum. Al- 
though using only a few bands in general does not yield 
accurate continuum spectra, it is sufficient for the con- 
vergence of exciton binding energies. Binding energies of 
Frcnkel excitons converge quickly with respect to the k- 
grid, while Wannier excitons require larger k-grids than 
usually seen in the literature. Although excitons are con- 
ventionally described as bound electron-hole pairs, only 
taking k-points near the F-point does not give a good de- 
scription for excitons in TDDFT, suggesting very strong 
collective character when represented with Kohn-Sham 



excitations. 

We test our formalism with several xc kernels. The 
LRC empirical formula, whose empirical parameter has 
been designed for reproducing the continuum spectrum, 
usually produces bound excitons as well, though the 
binding energies are generally too small. This is of course 
hardly surprising, because it is difficult to imagine how 
a single parameter could be sufficient to fit all aspects of 
the optical response. If one is interested in bound exci- 
tons rather than the continuum spectrum, the strength 
of the LRC kernel has to be increased. 

The bootstrap kernel is generally accurate for Frenkel 
excitons, while it produces Wannier excitons that are 
somewhat too weakly bound. On the other hand, the 
PGG kernel does not yield any bound excitons at all. 
Thus, at present we do not know of any simple, nonem- 
pirical xc kernel that produces accurate bound Wannier 
excitons in solids, xc kernels derived from many-body 
thcory^'^^'^'''^® may be expected to perform better than 
the kernels we have studied here, but they are signifi- 
cantly more complex. 

We also extended our formalism to triplet excitons, 
and derived a formula to generalize the bootstrap kernel 
to spin-dependent systems, with the hybrid kernel as a 
special case. This hybrid kernel yields the correct order of 
magnitude for singlet-triplet exchange splitting in some 
cases, but is in general not very quantitatively accurate. 
Thus, we have given a proof of principle that TDDFT 
is capable of producing reasonable exchange splittings in 
solids; the search for more accurate spin-dependent xc 
kernels for triplet excitons remains an important task for 
future investigations. 

In summary, we have shown that TDDFT shows con- 
siderable promise for treating excitonic effects, but more 
accurate multipurpose xc kernels for solids are needed, 
particularly for spin-dependent phenomena. Our ap- 
proach for directly calculating exciton binding energies 
will be convenient for facilitating such future develop- 
ments. 
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